% By Martin M. Andreasen, April 22 2010
% This function simulates the model when solved up to third order. The purning scheme is used.
% These codes are faster for a third order approximation because we use the kronecker
% product to eleminate loops
function [Y_sim,yields_sim,ER_sim,X_sim] = sim3rd(model,shocks,orderApp)

% Unfold model
gx  = model.gx;
%gxx = model.gxx;
%gxxx= model.gxxx;
gss = model.gss;
gssx= model.gssx;
GGxxtil  = model.GGxxtil;
GGxxxtil = model.GGxxxtil;

hx  = model.hx;
%hxx = model.hxx;
%hxxx= model.hxxx;
hss = model.hss;
hssx= model.hssx;
HHxxtil  = model.HHxxtil;
HHxxxtil = model.HHxxxtil;

eta = model.eta;
sig = 1;

numSim = size(shocks,2);    %Number of simulations

% Allocating memory
nMat       = length(model.by0);
Y_sim      = zeros(model.ny,numSim);
yields_sim = zeros(nMat,numSim);
ER_sim     = zeros(length(model.matConShortRate),numSim);
X_sim      = zeros(model.nx,numSim);

% The simulation
vecP0   = (eye(model.nx^2)-kron(model.hx,model.hx))\reshape((model.eta*model.eta'),model.nx^2,1);
xt      = (eye(model.nx)-model.hx)\(model.HHxxtil*vecP0+0.5*model.hss);
%xt         = zeros(model.nx,1);    

if orderApp == 1
   for t=1:numSim
      X_sim(:,t)      = xt;
      Y_sim(:,t)      = gx*xt;
      yields_sim(:,t) = model.byx*xt;
      ER_sim(:,t)     = model.rx*xt;
      
      % The states next period
      xt_p  = hx*xt + eta*shocks(:,t);
    
      % Updating xt
      xt  = xt_p;
   end
elseif orderApp == 2
   for t=1:numSim
      X_sim(:,t) = xt;
      AA = kron3(xt,xt);
      Y_sim(:,t)      = gx*xt        + GGxxtil*AA        + 0.5*gss;     
      yields_sim(:,t) = model.byx*xt + 0.5*model.byxx*AA + 0.5*model.byss;
      ER_sim(:,t)     = model.rx*xt  + 0.5*model.rxx*AA  + 0.5*model.rss; 
      
      % The states next period
      xt_p  = hx*xt + HHxxtil*AA + 0.5*hss + eta*shocks(:,t);
    
      % Updating xt
      xt  = xt_p;
   end
elseif orderApp == 3
   for t=1:numSim
      X_sim(:,t)      = xt;
      AA              = kron3(xt,xt);
      AAA             = kron3(xt,AA);
      Y_sim(:,t)      = gx*xt + GGxxtil*AA + GGxxxtil*AAA + 0.5*gss + 3/6*gssx*xt;
      yields_sim(:,t) = model.byx*xt + 0.5*model.byxx*AA + 1/6*model.byxxx*AAA + ...
                        0.5*model.byss + 3/6*model.byssx*xt;
      ER_sim(:,t)     = model.rx*xt + 0.5*model.rxx*AA  + 1/6*model.rxxx*AAA + ...
                        0.5*model.rss + 3/6*model.rssx*xt; 
               
      % The states next period
      xt_p  = hx*xt + HHxxtil*AA + HHxxxtil*AAA + 0.5*hss + 3/6*hssx*xt + eta*shocks(:,t);
        
      % Updating xt
      xt  = xt_p;
   end
end

end